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Abstract. A boundary-element alternating method, denoted herein as BEAM, 
is presented for two-dimensional fracture problems. This is an iterative method 
which alternates between two solutions. An analytical solution for arbitrary poly- 
nomial normal and tangential pressure distributions applied to the crack faces of 
an embedded crack in an infinite plate is used as the fundamental solution in the 
alternating method. A boundary-element method for an uncracked finite plate 
is the second solution. For problems of edge cracks a technique of utilizing “fi- 
nite elements” with BEAM is presented to overcome the inherent singularity in 
boundary element stress calculation near the boundaries. Several computational 
aspects that make the algorithm efficient are presented. Finally the BEAM is ap- 
plied to a variety of two-dimensional crack problems with different configurations 
and loadings to assess the validity of the method. The method gives accurate 
stress-intensity factors with minimal computing effort. 


1 Introduction 

Stress-intensity factors are fundamental fracture parameters that are needed to 
design structures against fatigue and fracture failures. In two-dimensional analysis, 
several methods are available in the literature to calculate the stress-intensity 
factors of cracked bodies. Several stress-intensity factor compendia [1-4] are also 
available. Recent research [5-13] revealed the potential of the alternating method 
to obtain stress-intensity factors in cracked bodies for which solutions are not now 
readily available. The alternating method developed and employed here is based 
on an earlier method known as the the Shwartz- Neumann alternating method 
[5]. The alternating method is an iterative numerical technique that alternates 
between two solutions to satisfy the required boundary conditions of the problem. 
One solution is a fundamental and is usually a continuum solution for a cracked 
infinite plate or solid. The second solution is provided by a numerical analysis 
such as finite- elements or boundary- elements of the uncracked body subjected to 



the same loading conditions. The method alternates between these two solutions 
to satisfy the required boundary conditions of the original problem: ~ 

The literature contains several papers on the alternating method which use 
the finite element method to obtain the second solution [6-12]. Recently the bound- 
ary element method was used in the alternating method instead of finite-element 
method [13]. The boundary element method (BEM) is attractive because with this 
method only the boundaries of the problem need to be modeled and hence, the 
modeling effort is considerably reduced. The purpose of this paper is to present 
such a boundary- element alternating method (BEAM} for two-dimensional mixed- 
mode crack problems. The method is thoroughly discussed and several attractive 
computational features of this method are highlighted. A procedure to combine 
“finite elements” with BEAM to overcome inherent singularities in stress calcula- 
tion with BEM is discussed. First, the basic analytical solutions for a crack in an 
infinite plate subjected to arbitrary normal and tangential pressure distributions 
over the crack faces are summarized. Second, a brief summary of the boundary- 
element method is presented to facilitate the presentation of the BEAM. Next, 
various computational aspects involved in the BEAM are discussed. Finally, sev- 
eral numerical examples are presented to illustrate the versatility of this method 
to obtain accurate mixed-mode stress-intensity factors for several cracked config- 
urations. 


2 Analytical solutions 


Use of the alternating method requires analytical solutions for an embedded crack 
in an infinite plate subjected to tractions on the crack faces. Consider then an 
infinite plate with a crack of length 2a as shown in Fig. 1. It is assumed that 
the tractions on the crack faces can be represented using an arbitrary polyno- 
mial pressure distribution, p yj applied normal to the crack faces and an arbitrary 
polynomial pressure distribution, applied tangential to the crack faces. The 
normal, ( p y ), and tangential, (p x ), pressure distributions applied over the crack 
faces take the polynomial form, 
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Figure 1 . Crack in an Infinite plate subjected to normal 
and tangential forces 
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where jV is the maximum degree of the polynomial functions and A n and B n are 
arbitrary constants. To obtain the stresses everywhere in the infinite plate due 
to the above arbitrary pressure distributions, the analytical solutions for a typical 
polynomial term ( x/a) n in Eq. (1) is obtained. The Westergaard stress function 
for any integer power n can be written as [10,11] 


V-(z) = < 


(z/a) 71 -| for odd n 

I ~( z / a ) n + -, — , S1 for even n 
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where z is a complex variable, z — x + iy and i = yf— 1 
The functions G^z) and G^z) are given in reference 10 as 
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for odd values of n with n — 2m + 1 


2 I 


*=0 


(4) 


and for even values of n with n = 2m. In Eqs. (3) and (4) the constants Ck are 
/ 1 for k = 0 

* “ i (iKt - l)(i - 2) - - - (i - fc + !) for * = 1,2,3... U 

The stress function of Eq. (2) is valid for both normal and tangential distributions 
of Eq. (1). 

The stresses at any general location z in the plate due to normal, p y , pressure 
distribution on the crack faces are [1] 


(*■*),, = SR(^) - yS(tf') 

= + ( 6 ) 

Similarly, the stress at any location z in the plate due to tangential, p z , pressure 
distribution, on the crack faces are 
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In equations (6) and (7), 9?( ) and £y( ) denote the real and imaginary parts, 
respectively, of the function in the parenthesis and 


,, _ # 
^ dz 


( 8 ) 


The stress-intensity factor weights, ( k w ) n , for mode / and mode II for a typical 
term (i/o) n in Eq. (1) can be obtained from the stress function as 


k w = v / 27r lim {\/z — a V’( z )} 
for the crack tip at x = a and 


k w — \Z2 k lim {(— i)y/z + a if>(z)\ 

z—> — a ^ } 


(9a) 


(96) 


for crack tip at x = — a. 

The stress-intensity factors for each of the polynomial terms ( x/a) n were 
computed and are presented in Table 1. 


3 The boundary element method 

In this section, the boundary-element method (BEM) is briefly outlined to facil- 
itate the presentation of the BEAM. In the absence of body forces the integral 
representation of the displacement at any internal point 1 j* for a two-dimensional 
domain can be written as [14] 
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Table 1 Stress-intensity factor weights for a crack in an infinite plate subjected to 
arbitrary normal (p v ) and tangential pressure (p x ) distributions of the form (x/ a) 


n 

(&t i>)n 

| 

1 

0 

i 

i 

1 

±1/2 

i 

2 

1/2 

i 

3 

±3/8 

1 

4 

3/8 


5 

±5/16 


6 

5/16 


7 

±35/128 


8 

35/128 



The postive and negative signs in this table refer to the crack tips 
at x = ±a, respectively. The negative values are meaningful only 
in the presence of additional forces which prevent crack closure. 


u l(Zj) “1“ J P*k( x jy£j) u k( x j) dT — J Pk( x j) dT (10) 

where u* k , p* k are components of displacements and tractions obtained from 
the fundamental solution of a unit point load in an infinite domain. Specifically, 
u *k( x j>€j) ls the displacement along the x/— direction at point x^, due to a unit 
point load along the x* — axis at point in an infinite domain of elastic linear 
isotropic material. Similarly, P*fc( x j^j) is the traction along the x/— axis on a 
line oriented at point Xj due to a unit point load along the x axis at point in 
an infinite domain. The quantity u *.(£., ) are the displacements at point £ j and so 
on. Equation (10) is the well-known Sormgliana’s identity and it gives the values 
of the displacements at any internal point in terms of the boundary values Uk,Pk 
[14]. 

In the BEM the boundary (or boundaries) of the domain is discretized into 
N e number of boundary elements. The displacements Uk and tractions pk within 
each of these elements are assumed to vary in a constant or linear or quadratic 
manner. The discretized form of Eq. (10) is written as 
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A boundary integral expression can be obtained by taking Eq. (11) to the bound- 
ary. This equation when applied to different points on the boundary produces a 
system of equations that are of the form 


{H]{U} = [<5]{r} (12) 

where {£/} and {T} are the displacements and tractions, respectively, at all the 
nodes on the boundary I\ For a mixed boundary value problem, the displace- 
ments {U p } are prescribed on some parts of the boundary and tractions { T p } are 
prescribed on the remaining parts of the boundary. Note that at any point on 
the boundary either displacements or tractions are known. Equation (12) can be 
rearranged to group all the unknowns on the left hand side (LHS) and all the 
known quantities on the right hand side (RHS). This rearrangement leads to 


mm = ig]{f} 


(13) 


where {A} is the column vector of all the unknown displacements and tractions 
and {F 7 } contains the values of the known displacements and prescribed tractions. 
Note that unlike finite-element method, the matrix [H] in Eq. (13) is unsymmetric 
and is fully populated. The unknown displacements and tractions can be obtained 
solving Eq. (13) as 


{X} = [H]-' [G]{F} 


(14) 


Once the nodal displacements and tractions are known, the displacements, {u}, 
and stresses, {<r}, at any internal point can be calculated using the Somigliana’s 
identity as, 
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W = £ m j {T e } j - £ [D 2 \j{U e }j 

j=l j = 1 

N« N t 

M = E l s >]y{ r '}y - E IH^}; 

j -1 1=1 

where iV c is the number of boundary elements and {T c } and { U e } are nodal trac- 
tions and displacements of the j th element. The matrices [Di], [£> 2 ]? [Si] and [ 52 ] in 
Eqs. (15) are obtained from integrating the integrals involved in displacement and 
stress calculations [14]. However, when the internal point approaches the bound- 
ary the integrals involved in the stresses and displacements at the internal point 
become hyper-singular [14] and hence give rise to numerical errors. A procedure 
that circumvents this difficulty when the internal point approaches the boundary 
and obtains accurate stresses and displacements is presented later in this paper. 


4 Boundary element alternating method 

The procedure used in BEAM is very similar to that used in the finite-element 
alternating method in references [6-12] The procedure is illustrated in Fig. 2, and 
is explained below. 

Step 1 : Analyze the same configuration and loading as the given problem, but 
without the crack using the boundary element method. The boundary element 
solution gives the displacements and tractions of all the nodes on the boundaries 
of the uncracked body. (Note that the line y = 0 corresponds to the line on which 
the crack is present in the original problem.) 
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Figure 2. Procedure used in the alternating method 


Step 2 : Obtain the stresses, a y and cr zy , on the line y = 0 for —a < x < a 
using Eq. 13. These stresses determine the normal and tangential tractions p y and 
p z , respectively on this line y = 0. 

Step 3 : If both the normal and tangential tractions are negligibly small (i.e. 
smaller than a prescribed tolerance level) stop the algorithm and calculate the 
sum of the stress-intensity factors computed so far. If either normal or tangential 
or both tractions are not negligible go to Step 4. 

Step 4 : For the crack faces to be traction free, the crack face normal tractions 
(p y ) and tangential tractions (p z ), computed in Step 3 must be removed. (In 
Fig. 2, for clarity, the tangential tractions on the crack faces are not shown.) 
Therefore, the negative of crack-face normal tractions (i.e. p y = — p^) and crack 
face tangential tractions (i.e. p x = —p x ) are applied to the analytical solution. 
The tractions p ^ and p ^ are now expressed in polynomial form as 

p? = E *.(*/•)" = mV) 

n=0 

(16) 

p? = E n n (x/*r = {p) t {b} 

ti = 0 
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where 


{P} T = { 1 , (*/“)» ( a: / 0 ) 2 > (*/ a )™} 

The coefficients {A} and {B} in Eq. (16) are calculated using a least square 
procedure [6] by 


{A} = 

{B} = 

where 

> £i =/_ 

w -s: 


[■ e r [d] 

[E]-' [C] 

(17) 

{P}{P} T dx 

(18) 

{P}Pyi X ) dx 

(19) 

{P}Px( X ) dx 

(20) 


The integrals in Eqs. (18)-(20) can be computed easily by using numerical inte- 
gration, such as Gaussian quadrature, because the discrete numerical values of 
p^(x) andp^(x) can be calculated at Gaussian points from the boundary element 
solution. 

Step 5 : Once the coefficients {^4} and {5} are determined from Eq. (17), 
equations (9) are used to calculate the stress-intensity factors for the current j th 
iteration as 


xi-E 

n=° (21) 

N v ’ 

Kjl = (ku>) n -®n 

n=0 

where (k w ) n are the stress-intensity factor weights given in Table 1 for each of the 
polynomial functions. 
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Step 6 : The crack-face normal p^ and tangential p * tractions obtained in 
Step 4 create tractions on all the boundaries of the region of interest. The stresses 
at any point z = x + iy , on the boundary can be obtained using Eqs. (2)-(8) for 
each of the polynomial functions as 


{< 7 } = [M]{A} + \L]{B} (22) 

where 

M T = {°z a v <r*y} ( 23 ) 

In Eq. (22), [M] and [L] are the stresses at any point z on the boundary due to 
unit values of each of the polynomial normal and tangential pressure distributions, 
respectively. Note that the matrices [M] and [L] are available from Eqs. (6) and 

( 7 ). 

The tractions, T x and T y , in the x— and y— directions, respectively, can be 
calculated using the stresses on the boundary as, 


m = mm 


(24) 


where 


m T ={T, Ty}, 


Tlx Tly 

0 n y n x 


(25) 


In Eq. (25), n x and n y are the direction cosines of the normal to the boundary 
with respect to x— and y— axes, respectively. 

Step 7 : To satisfy the traction-free boundary conditions on the external 
boundaries, the tractions created on the these boundaries due to the residual 
pressures on the crack face in Step 6 need to be removed. Therefore, the negative 
of these tractions calculated at the nodes in Step 6 are applied to the boundary 
element model of the uncracked plate. The tractions calculated in Eq. (22) at each 
node of the model are assembled to form a global vector as 


11 



(26) 


{T} = -[ GMA] - [G,]{fl} 

where [G n ] and [Gt] are the assembled matrices obtained using Eqs. (22), (24) and 
(25) for each node in the model. 

The nodal tractions in Eq. (26) along with the original displacement boundary 
conditions are prescribed as tractions and displacements to the uncracked body. 
This is the start of the next iteration. This iterative procedure is continued until 
the crack face tractions in Step 3 are negligibly small. In the converged solution, 
the mode I and mode II stress-intensity factors are simply the sum of the stress- 
intensity factors from all the iterations. The convergence criteria were formulated 
in terms of the integral of the residual pressures on the crack faces and are de- 
scribed in reference 10. Here, similar convergence criteria are used for both normal 
and tangential tractions simultaneously. 


5 Some Computational Aspects 

In the boundary element part of the method parabolic straight elements with 
3-nodes were used to model the uncracked boundaries of the body. Several com- 
putational aspects of the alternating method which are generic to any numerical 
method used to analyze the uncracked body are presented in [6,7,10,12]. Other 
computational aspects pertinent to the BEAM are summarized below. 

5.1 Residual Pressure Calculation 

The residual pressures fit over the crack center line (where the crack is present 
in the original problem) and each iteration requires evaluation of the integrals in 
Eqs. (18)-(20). To evaluate these integrals the values of stresses at quadrature 
points are required. These stresses at quadrature points are calculated by using 
Eq. (15) and assembling contributions from all the boundary elements in the struc- 
ture. This needs to be performed for each iteration since, the displacements and 
tractions used in Eq. (15) differ from iteration to iteration. However, considerable 
computational efficiency can be achieved by a technique described below. 

For any given crack length a , the coordinates of the quadrature points used 
in Eqs. (18)-(20) do not vary from iteration to iteration. Hence, before starting 
the iterative process the matrices [Si] and [S 2 ] in Eq. (15) for each of the bound- 
ary elements are calculated at the quadrature points. These matrices are then 
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assembled to form the global matrices [S* g ] and [S^] for each of the quadrature 
points. Then the stresses at a quadrature point can be evaluated from, 


w = [s./im* - (27) 

where {U} k and {T} fc are the global nodal displacement and traction vectors for 
the k th iteration. The stresses at all the quadrature points for each iteration are 
obtained by simple multiplication involved in Eq. (27). This procedure, however, 
requires storage of the matrices [S* g ] and [5^] at each of the quadrature points. 


5.2 Calculation Of Nodal Tractions And Displacements 

The unknown nodal tractions and displacements for any k th iteration are obtained 
by combining Eqs. (13) and (26) as 

{X} t = [<?]{A}, + [*]{£}, (28) 


where 


[Q] = -[IT 1 [G][<?n] 

[R] = -[H]- 1 [G\[G t ] 


(29) 


and [G n ] a ^d [Gt] are formed from the matrices [G n ] and [G t ] of Eq. (26) and the 
prescribed displacement boundary conditions. 

As seen from Eq. (28), the unknowns {X}*, for the k th iteration depend only on 
the polynomial coefficients {A}* and and the matrices [Q] and [R]. The 

matrices [H] } [£?], [C? n ], [Gt] and hence [Q] and [R] do not vary in the iterative 
process. Therefore, the assembled LHS matrix [H] needs to inverted only once 
and the matrices [Q] and [i2] need to be computed only once using Eq. (29). 
This can be performed conveniently before the start of the iterative process. For 
each iteration, after the residual pressure fit is made to evaluate {A}*, and {5}*, 
the unknowns {X}* are immediately available from Eq. (28) by simple matrix 
multiplications. 
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Figure 3. "Finite elements" near the edge for accurate stresses 


5.3 Stresses and Displacements Near Boundary Points 

For edge crack problems, the residual polynomial pressures over the crack line 
requires estimation of stresses from the BEM portion of the analysis at points very 
near the boundary where the crack emanates. As mentioned earlier, as the internal 
point approaches the boundary, the integrals involved in Eq. (15) for displacements 
and stresses at internal points become hyper-singular [14] and because of this 
singularity the stresses and displacements calculated very near the boundary, in 
general, are inaccurate. This can be overcome by integrating the integrals in 
Eq. (15) in closed form. However, this is not feasible for higher order boundary 
elements. Therefore, when the internal points are very close to the boundaries 
special procedures need to adopted to obtain accurate stresses. One such procedure 
is described below. 

The displacements calculated near the boundary points are more accurate 
than the stresses calculated at the same point. This is because, the singularity in 
Eq. (15) for displacement calculations near the boundary is one order less than the 
singularity in stress calculation. This fact is utilized in an extrpolation procedure. 
In this procedure two “finite elements” are constructed as shown in Fig. 3, near the 
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point 63 on the edge from which the crack emanates. In Fig. 3 , the crack is along 
the line 63 x. The nodes 61,62,63,64,65 are the nodes in the boundary element 
model and j and ( j + 1) are the elements near the “crack line” of the uncracked 
body. The two “finite elements” (Fi and F 2 in Fig. 3 ) are constructed such that 
they (a) are like 8-node parabolic finite elements, (b) axe physically located on 
either side of the crack line and share a common side (637i 2 7i 6 )- the crack line, (c) 
share nodes that belong to the boundary element idealization of the uncracked 
body (61,62,63 for the element F 1 and 63,64,65 for element F 2 ) and (d) penetrate 
a distance A into the uncracked body. 

The displacement at the nodes rii,n 2 , . . . ,7ig in the two “finite elements” can 
be computed by using Somigliana’s identity. The displacements at the bound- 
ary 61,62, ... ,65 are known from the boundary element solution of the uncracked 
body. Thus the displacements at all the 8-nodes of each of the finite elements F 1 
and F 2 are known. Following well known isoparametric finite element procedures 
the stresses at the 2 x 2 Gaussian points are calculated in each of the finite ele- 
ments. The 2x2 Gaussian stresses are then extrapolated to the crack center line 
[ 10 , 15 , 16 ]. In general, the crack line stresses thus obtained from the top element 
F\ are not identical to those obtained from the bottom element F 2 . Therefore, the 
extrapolated stresses from Fi and F 2 are averaged and are used in the residual fit. 

For points on the crack line whose distances from the edge are greater than 
A (like points mi,m 2 , . . .7724 in Fig. 3 ), the stresses can be calculted using the 
Somigliana’s identity , since the boundary element stresses are accurate at these 
points. 

In this procedure, the critical parameter is distance A, the length of the 
“finite elements”. Numerical experimentation suggests that the length A should 
be equal to or up to 1.5 times the length of the smallest element near the edge 
in the boundary element model of the uncracked body. This procedure appears 
to demand extra computations. However, the extra computing is well worth the 
effort because the stresses computed by this procedure are much more accurate 
than those obtained by conventional boundary element procedures. 

5.4 Consistent Tractions due to Residual Pressure Distribution 

In step 6 of the BEAM, the tractions on all the boundaries due to the residual 
pressure distributions p x and p y on the crack faces are determined. If the variation 
of the tractions over an element is higher than a quadratic, then directly using 
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the nodal values of the tractions computed from Eq. (24) in the boundary element 
equations will be inconsistent. This is because the traction computed at any point, 
I, on the boundary element j (see figure 4) using Eq. (24), T(s) is not the same 
as the traction computed T/ at the same point by using 


Tj — $\T\ + $2^2 + ^3^3 (30) 

where Tj , T 2 , T 3 are the nodal tractions calculated by Eq. (24) and $i,4>2,$3 are 
the quadratic shape functions of element j. The error E(s), at any point s on the 
j th boundary element is 


E(s) = Ti — T(s) 


(31) 


where {$} J = { $2 *3 }J and {T}J = { Tj T 2 T 3 }J 

The tractions {T} ■ that minimize the error E(s ) can be calculated by using 
a least square procedure as 


Element Q+1) 



Element (j-1) 


Figure 4. Tractions on element j on the boundary due to 
crack face pressure loading 
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(’’ «{4>} T <fcW = f’ {4.} t T(s)ds 
0 JJo 


( 33 ) 


where lj is the length of the m j th element. The Eq. (33) can be concisely written 
as 


l i \ 

{$} T(s)ds ) (34) 

The integrals involved in Eq. (34) can be carried out numerically because discrete 
values of T(s) on the elements are available. The matrix [W] -1 can be written 
down explicitly for a 3-node quadratic boundary element as 






-3/2 

9/4 

-3/2 


3 

-3/2 

9 


(35) 


5.5 Edge Cracks 

The analytical solution used in BEAM is for an embedded crack in an infinite 
plate having two crack tips. But in edge cracks problems, only one crack tip 
exists. Hence to use BEAM for edge crack problems, a fictitious crack tip needs 
to be defined. The fictitious crack tip is usually assumed to be at x = —a. While 
for 0 < x < a, the residual pressures are computed from the boundary element 
solution, the residual pressures on the fictitious part of the crack —a < x < 0 need 
to be defined. In this paper, a simple constant pressure distribution was used over 
the fictitious part as discussed in reference 12. 




(a) Without shift 


► X, X' 


Figure 5. Orgln shifting to account for the position of the fictious crack tip 
for edge cracks from internal boundaries 

If the crack length is large, the fictitious crack tip can penetrate the bound- 
ary opposite the point of origination of the actual crack. One such example of 
a long edge crack from a circular hole is illustrated in Fig. 5(a). There are two 
difficulties associated with long edge cracks and penetration of the fictitious part 
of the crack. First, pressure fitting needs to be carried out on the part of the 
fictitious crack that is in the penetrated region. This requires complicated bound- 
ary element idealization and the construction of “finite elements” at two places 
on the boundary. The second difficulty pertains to the use of the same boundary 
element model for all crack lengths. Although only the uncracked body is mod- 
eled, the boundary element idealization is usually made finer in the region of the 
crack and coarse elsewhere. For example, region AB in Fig. 5(a) is usually a fine 
mesh region and BC A is a coarse mesh region. When the cracks are short the 
fictitious crack tip, the one at (— a,0), is closer to the fine mesh region, boundary 
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AB in Fig. 5(a). In contrast, when the crack is long, the fictitious crack tip is near 
the region where the mesh is coarse, like boundary BC . Because of the square 
root singularity at the crack tips, for larger crack lengths the fictitious crack tip 
creates large tractions on the elements in the coarse mesh region. The iterative 
process may not be able to erase these tractions created in the elements of the 
coarse region. In such situations, the iterative process either oscillates or diverges. 
The current algorithm can be modified to make it more robust and allow the same 
boundary element model to be used for all crack lengths. The modification is to 
shift the origin such that the fictitious crack tip is positioned in such a way that 
its influence is felt most in the region of the fine mesh. Figure 5(b) shows such 
an origin shift by an amount ao so that the fictitious crack tip is close to the fine 
mesh region. This method is termed herein, the origin-shifting method. Details 
of the origin shifting method can also be found in reference [12]. 

These modifications to the BEAM facilitate the analysis of several crack 
lengths in a single computer run and require only one inversion of the system 
matrix in Eq. (14). Furthermore, two different coordinate systems are conveneint 
for the uncracked body analysis and the analytical solution. A global coordinate 
system is usually convenient for the uncracked body analysis, while a local (crack) 
coordinate system is convenient for the analytical part of the method. Appendix- 
A present the transformations that are necessary to be able use two different 
coordinate systems. 


6 Results and discussion 

To evaluate the effectiveness of the BEAM in mixed-mode problems, the method is 
applied to several cracked configurations and loadings for which accurate solutions 
are available in the literature. Both embedded and edge cracks in plates with and 
without stress concentrations (due to notches and holes) were analyzed. The 
normalized stress intensity factors are compared with those from the literature [4] 
for each of these problems analyzed. The percent difference in the tables presented 
in this paper is defined as 


percent difference — 


(BEAM value-Reference value| 
Reference value 


In all the problems analyzed, a plane strain state was assumed with a Poisson’s 
ratio of 0.3. 
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Experience with FEAM [12] showed that the degree of polynomial of value 
N = 5 yields accurate stress-intensity factors. Hence N = 5 was used for all the 
problems presented in this paper. Several crack lengths were analyzed with a single 
boundary element model for the uncracked body in a single computer run. All the 
boundary element idealizations were performed using straight 3-noded quadratic 
elements. When the prescribed tractions at a boundary point are discontinuous 
(such as a corner), double nodes were used near that point. 



(a) Remote uniform tension 
loading 


(b) Remote parabolic loading 


Figure 6: Embedded slant crack problems analyzed 
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6.1 Embedded Slant Crack Plate 


Remote Uniform Tensile Loading — An embedded slant crack specimen with re- 
mote uniform tensile stress is shown in Fig. 6(a). The plate was modeled with 10 
quadratic boundary elements and had 24 nodes. The sides X = were modeled 
with four equal length elements and the sides Y = izH were modeled with one 
element. Two crack inclination angles ( 6 = 45° and 75°) were considered. The 
boundary element models used are shown in Fig. 7. All the external boundaries, 
AB^BC, CD and DA y in Fig. 7 were made stress free in the iterative part of the 
algorithm. Eight crack lengths with ( a/W ) ratios varying from 0.1 to 0.8 were 
analyzed in a single computer run. Stress-intensity factors were calculated for all 
crack lengths at both crack tips ( x = ±a). The stress-intensity factors at both 
the crack tips, as expected, were computed to be identical (to machine accuracy). 
Tables 2 and 3 present the normalized stress-intensity factors obtained for 6 = 45 ° 
and 75° respectively. Reference results from the literature are included in these ta- 
bles for comparison. Excellent agreement is obtained for all crack lengths between 
the two sets of results even with a very coaxse boundary element idealization. The 
maximum difference is about two percent for both crack inclination angles. Three 
to five iterations were needed to achieve converged solutions. 



Figure 7. Boundary element model for slant embedded crack problems 
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Table 2 Comparison of normalized stress-intensity factors, F j and Fu, for embed- 
ded slant cracks in a rectangular plate subjected to remote uniform tensile stress. 
(Fig. 6(a ),(H/W) = 2.0, 9 = 45°) 


(a/WO 

Fi 

= x,/[ 

Sy/iTa ] 

Fu 

= K n / [ S ^ ] | 

BEAM 

value 

Reference Percent 

value difference 

BEAM 

value 

Reference 

value 

Percent 

difference 

0.1 

0.505 

0.505 

0.0 

0.502 

0.502 

0.0 

0.2 

0.518 

0.518 

0.0 

0.507 

0.507 

0.0 

0.3 

0.540 

0.541 

0.05 

0.515 

0.516 

0.39 

0.4 

0.571 

0.572 

0.18 

0.527 

0.529 

0.57 : 

0.5 

0.609 

0.613 

0.16 

0.543 

0.546 

0.73 

0.6 

0.658 

0.661 

0.45 

0.563 

0.567 

0.71 l 

0.7 

0.717 

0.721 

0.14 

0.589 

0.595 

1.00 : 

0.8 

0.785 

0.795 

0.25 

0.623 

0.630 

l.i : 


Table 3 Comparison of normalized stress-intensity factors, F j and Fu, for embed- 
ded slant cracks in a rectangular plate subjected to remote uniform tensile stress. 

(Fig. 6(a), (H/W) = 2.0, 9 = 75°) 


1 

i 

Fi = Kij [ Sy/iTa ]. 

Fu = Ki/ [ Sy/iTa ] 


a 



— 


= 

(a/W) 

BEAM 

value 

Reference 

value 

Percent 

difference 

BEAM 

value 

Reference 

value 

Percent 

difference 

- 

3 

0.1 

0.068 

0.068 

0.0 

0.252 

0.252 

0.0 


m 

0.2 

0.070 

0.070 

0.0 

0.256 

0.256 

0.0 


5 

0.3 

0.074 

0.074 

0.0 

0.263 

0.263 

0.0 



0.4 

0.078 

0.078 

0.0 

0.272 

0.272 

0.0 


2 

0.5 

0.083 

0.084 

1.2 

0.282 

0.283 

0.35 

— 

3 

0.6 

0.089 

0.090 

1.1 

0.293 

0.294 

0.34 

- 

i__ 

0.7 

0.094 

0.096 

2.1 

0.305 

0.306 

0.33 

H 

m 

0.8 

0.101 

0.102 

0.98 

0.317 

0.319 

0.32 

— 
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Table 4 Comparison of normalized stress-intensity factors, Fj and Fjj, for an em- 
bedded slant crack in a rectangular plate subjected to remote parabolic tensile 
stress. (Fig. 6(b), (H/W) = 2.0, 0 = 45°) 


(a/W) 

Ft 

= Ki/ [ Sy/ifa ] 

Fu 

= Kul [ Ss/ttcl ] 

BEAM 

value 

Reference 

value 

Percent 

difference 

BEAM 

value 

Reference 

value 

Percent 

difference 

0.25 

0.345 

0.344 

0.29 

0.361 

0.363 

0.55 

0.20 

0.338 

- 

- 

0.359 

- 

- 

0.30 

0.353 

- 

- 

0.365 

- 

- 

0.40 

0.375 

- 

- 

0.373 

- 

- 

0.50 

0.403 

- 

- 

0.384 

- 

- 

0.60 

0.438 

- 

- 

0.398 

- 

- 

0.70 

0.479 

- 

- 

0.415 

- 

- 

0.80 

0.528 

- 

- 

0.436 

- 

- 


Remote Parabolic Tensile Loading — Instead of remote uniform tensile loading, a 
parabolic tensile loading was applied to the slant crack specimens with a crack 
inclination angle of 45° as shown in Fig. 6(b). The applied stresses S(£) at a point 
on the top and bottom boundaries of the plate were calculated using the equation 

s({) = s[i.o-(^) 2 ] 

where —W < t > W . The boundary element model used for this problem was 
same as the previous problem (see Fig. 7). The mode I and mode II stress- 
intensity factors are presented in Table 4 along with the reference value. Reference 
results are available only for a crack length of ( a/W ) = 0.25 and again excellent 
agreement with the reference value was obtained with a very coarse boundary ele- 
ment idealization. A converged solution was obtained with three to five iterations. 


These examples demonstrate that the BEAM gives accurate mixed-mode 
stress-intensity factors for problems of embedded cracks in plates with rectilin- 
ear boundaries. The method also effectively handled various loading conditions. 


23 


In the next section the BEAM is applied to problems involving edge cracks in 
mixed-mode situations. 



(a) Slant edge crack- 
Remote tension 
loading 



(b) Slant edge crack- 
Remote bending 
loading 



(c) Slant edge crack from a 
circular hole 


Figure 8. Slant edge crack problems analysed 


6.2 Slant Edge Cracks in Plates 

Remote Uniform Tensile Loading — A slant edge cracked plate with remote uni- 
form uniaxial loading is shown in Fig. 8(a). A crack inclination angle of 45° 
was used. The boundary element model with 56 nodes and 26 elements used in 
the analysis is shown in Fig. 9. Two finite elements were constructed near the 
point from which the crack is emanating as shown in Fig. 9. The distance A as 
discussed earlier was selected to be about 1.5 times the length of the smallest 
boundary- element near the edge. Crack lengths varying from ( a/W ) = 0.3 to 
0.6 were considered. The stress-intensity factors obtained using this model are 
presented in Table 5 along with the reference solutions. The maximum difference 
between the results from this analysis and the reference solution is about 3 percent. 
About 10 to 12 iterations were needed to obtain converged solutions. 

Even though the shift of the origin is not necessary for the edge cracks ema- 
nating from external boundaries, the origin-shifting method is used in the above 
analysis, since numerical experimentation showed that the origin shift gave better 
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results than the one without origin shifting. This may be due to the effect of 
the fictitious crack tip is minimized in the origin shifting method. In the above 
analysis the fictitious crack tip was positioned at a distance ( x/W ) = —0.3 from 
the free edge. 



Figure 9. Boundary element model for edge crack problems 


Remote Bending Loading — A slant edge cracked plate with remote bending is 
shown in Fig. 8(b). The boundary element model is shown in Fig. 9. The bend- 
ing moment was simulated by the corresponding bending stress applied over the 
edges as shown in Fig. 8(b). Crack lengths ranging from ( a/W ) = 0.3 to 0.6 
were considered. The results are presented in Table 6 with reference results from 
the litreature. Here again excellent agreement was obtained with the reference 
solutions. 
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Table 5 Comparison of normalized stress-intensity factors, Fj and Fjj , for slant 
edge cracks in a rectangular plate subjected to remote uniform tensile stress with 


origin shift. 

(Fig. 8(a), 

= 1.0, (H 2 /W) = 1.5, e 

= 45°) 


| 

(«/w0 

Fi 

= Kjj [ Sy/rTa ] 

Fu 

= Kii/ [ Sy/ira ] 

i 

BEAM 

value 

Reference 

value 

Percent 

difference 

BEAM 

value 

Reference 

value 

Percent 

difference 

0.10 

0.728 

- 

- 

0.400 

- 

- 


0.20 

0.792 

- 

- 

0.425 

- 

- 

l 

0.30 

0.889 

0.879 

i.i 

0.465 

0.450 

3.3 

- 

0.35 

0.943 

0.940 

0.32 

0.490 

0.473 

3.6 

_ 

0.40 

1.009 

1.011 

0.20 

0.518 

0.505 

2.6 

= 

0.45 

1.080 

1.094 

1.3 

0.551 

0.536 

2.8 

7 

0.50 

1.184 

1.190 

0.50 

0.590 

0.574 

2.8 


0.55 

1.311 

1.301 

0.77 

0.636 

0.616 

3.3 

- 

0.60 

1.434 

1.437 

1.2 

0.699 

0.674 

3.7 

- 


Table 6 Comparison of normalized stress-intensity factors, Fi and Fji , for slant 
edge cracks in a rectangular plate subjected to remote bending stress. (Fig. 8(b), 
(Hi/W) = 1.0, (H 2 /W) = 1.5, 6 = 45°) 


Fj = Kij [ S \/Tr~a } F u = K n j [ ] 


9 

(a/W) 

BEAM 

Reference 

Percent 

BEAM 

Reference 

Percent 

i 

value 

value 

difference 

value 

value 

difference 

1 

| 

0.10 

0.674 

- 

- 

0.360 

- 

- 

i 

0.20 

0.673 

- 

- 

0.345 

- 

- 


0.30 

0.690 

0.685 

0.73 

0.334 

0.315 

6.0 

§ 

0.35 

0.705 

0.705 

0.0 

0.332 

0.319 

4.1 

i 

0.40 

0.728 

0.736 

1.1 

0.330 

0.320 

3.1 

I 

0.45 

0.756 

0.758 

0.26 

0.331 

0.324 

2.2 

1 

- 

0.50 

0.789 

0.795 

0.76 

0.340 

0.329 

3.3 

1 

0.55 

0.827 

0.832 

0.6 

0.341 

0.337 

1.2 

1 

0.60 

0.872 

0.889 

1.9 

0.352 

0.344 

2.3 
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6.3 Slant Edge Crack Emanating From A Hole In A Rectangular Plate 

An edge crack emanating from a hole with radius R in an infinite rectangular 
plate and loaded with remote uniform tension is shown in Fig. 8(c). The crack 
inclination angle 0 was selected as 30°. Since the reference solution available for 
this problem is of an infinite plate, a plate with dimensions ( H/W ) = 1.0 and 
(W/R) = 24 was considered. The outer boundaries of the plate were modeled 
with 12 elements as shown in Fig. 10 (a). The hole boundary was modeled with 
78 boundary elements as shown in Fig. 10(b). ( Note that the middle node of the 
boundary elements are not shown in Fig. 10(b).) Here again two finite elements 
were constructed near the crack mouth as shown in Fig. 10(c). As described in 
[12], for all crack lengths with (a/ R) > 0.25, the fictitious crack tip was positioned 
at ( x/R ) = —0.25. The crack lengths ranged from ( c/R ) = 0.1 to 1.0. (Note that 
c is the projected length of the crack along the horizontal, see Fig. 8(c).). The 
stress-intensity factors obtained with this method are presented in Table 7 along 
with the reference solution. Excellent agreement was obtained with the reference 
solutions. Only 10 to 14 iterations were needed to obtain converged solutions. 



(a) Modal for the plate 


(b) Model for the circular hole 


Figure 10. Boundary element model for an edge crack from a circular hole problem 
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Table 7 Comparison of normalized stress-intensity factors, Fj and F/j, for slant 
edge cracks emanating from a circular hole in a rectangular plate subjected to 
remote uniform tensile stress. (Fig. 8(c), ( H/W ) = 1.0, 0 = 30°) 


(c/R) 

Fi 

= if// [ Sy/iTc } 

Fn 

= Knl [ S^Z ] 

BEAM 

value 

Reference 

value 

Percent 

difference 

BEAM 

value 

Reference 

value 

Percent 

difference 

0.10 

2.501 

2.501 

0.0 

0.750 

0.701 

7.0 

0.15 

2.328 

(2.300)* 

1.2 

0.677 

(0.660) 

2.6 

0.20 

2.165 

2.140 

1.2 

0.617 

0.610 

1.2 

0.30 

1.866 

(1.880) 

0.75 

0.548 

(0.550) 

0.36 

0.40 

1.645 

(1.680) 

2.1 

0.519 

(0.520) 

0.19 

0.50 

1.501 

1.519 

1.2 

0.506 

0.511 

0.98 

0.75 

1.226 

(1.240) 

1.1 

0.490 

(0.500) 

2.0 

1.00 

1.060 

1.090 

2.8 

0.491 

0.508 

3.3 


* Values in parantheses are interpolated values from tables in reference 4 (pages 
259 and 261). 


7 Concluding remarks 

A boundary element alternating method (BEAM) for two-dimensional, 
mixed-mode fracture problems is presented. The analytical solution for arbitrary 
normal and tangential pressure distributions on the faces of the crack in an infinite 
plate is used as the fundamental solution. In the numerical part, the boundary 
element method is used to model and analyze the uncracked body. Details of 
the implementation of the algorithm, together with a variety of computational 
aspects of the method, were presented. An origin-shifting method is presented 
that is particularly useful in the analysis of edge cracks emanating from the in- 
ternal boundaries. This method of origin-shifting was also found to give better 
results even for edge cracks emanating from internal boundaries. Also for edge 
cracks problems, a procedure which utilizes “finite elements” and finite-element 
type interpolations and stress evaluation was developed to calculate the stresses 
accurately near a boundary. This procedure yielded stresses that are of better ac- 
curacy than the stresses calculated by the traditional boundary element method. 
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The BEAM was applied to several mixed-mode problems to evaluate its ef- 
fectiveness. These numerical examples showed that the BEAM requires very little 
modeling effort and yields accurate mixed-mode stress-intensity factors. Three 
to four iterations were necessary to yield accurate and converged stress-intensity 
factors, while edge crack configurations needed about 10-15 iterations. With this 
method, several crack lengths can be analyzed in a single computer run and hence 
it can be used to economically obtain stress-intensity factors over a range of crack- 
lengths. 
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Appendix - A 

Global-local Transformations in the BEAM 

This appendix describes the coordinate and stress transformations that are 
required when the global coordinate system does not coincide with the local (crack) 
coordinate system. 

The analytical part of the alternating method assumes that the crack is along 
the y = 0 line (i.e. x — axis). While the crack coordinate, (x — y), system is 
convenient for the analytical solution, it is cumbersome for the boundary element 
part of the method. For example in Fig. A-l the global, (X — F), coordinate system 
is obviously convenient for the BE part of the method and not the x — y crack 
coordinates system. These two coordinate systems require stress and coordinate 
transformations in the BEAM. These are discussed below. 


rTTT^ 



Figure A-1, Global-local coordinate systems In BEAM 


Before performing the Step 2 of the BEAM, the stresses, (Tx,^Y and < txy, 
obtained from the uncracked body solution need to transformed to the local (crack) 
coordinate system as, 
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(A-l) 


f 1 1 

[ 

II 

b 

cry 

{<Txy ) 1 

l axY , 


In Eq. (A-l), [A] is the transformation matrix given by 





sin 9 cos 9 


sin 2 9 
cos 2 9 
sin 9 cos 9 


sin 29 
— sin 29 
cos 29 


(A -2) 


where 6 is the crack inclination angle measured from the X — axis of the global 
coordinate system as shown in Fig. A-l. The transformed stresses <x x ,<x y and cr xy 
are then used to obtain the normal and tangential tractions on the crack faces. 

Step 6 of the alternating method, on the other hand, requires an inverse trans- 
formation. Recall that in Step 6 of the method the stresses on all the boundaries of 
the body due to crack face pressures p ^ and p ^ are calculated. From these stresses 
the tractions on the boundaries are calculated. First to calculate the stresses at a 
point (Xb, Y&) on the boundary, the coordinates of this point in the local system 
are calculated as 



— sin 0 1 j X b -X c \ 
cos 0 \{Y b -Y c j 


04 - 3 ) 


where X C ,Y C are the global coordinates of origin of the local coordinate system. 
Then the complex variable at this point can be formed as, z b = x b + iy b . This 
complex variable is then used in Eq. (22) to evaluate stresses in the local coordinate 
system. These stresses are then transformed to the global system using 


( tr x 


\ aY 

> = [^] _1 { ^ » 

{ &XY ) 

i l a xy j 


04 - 4 ) 


These global stresses are then used in Eq. (24) in Step 6. 

With the above definitions, however, Eq. (25) need to be redefined as, 
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M-5) 


{T} t = {T x TV} 


M = 


Tlx 

0 


Tly 0 
Tly Tlx 


where nx and ny are the direction cosines of the normal to the boundary with 
respect to X— and Y— axes, respectively. 
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